Primordial Non-Gaussianity and Reionization 
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The statistical properties of the primordial perturbations contain clues about the origins of those 
fluctuations. Although the Planck collaboration has recently obtained tight constraints on primor- 
dial non-gaussianity from cosmic microwave background measurements, it is still worthwhile to mine 
upcoming data sets in effort to place independent or competitive limits. The ionized bubbles that 
formed at redshift z ~ 6 — 20 during the Epoch of Reionization are seeded by primordial overden- 
sities, and so the statistics of the ionization field at high redshift are related to the statistics of the 
primordial field. Here we model the effect of primordial non-gaussianity on the reionization field. 
The epoch and duration of reionization are affected as are the sizes of the ionized bubbles, but these 
changes are degenerate with variations in the properties of the ionizing sources and the surrounding 
intergalactic medium. A more promising signature is the power spectrum of the spatial fiuctuations 
in the ionization field, which may be probed by upcoming 2i cm surveys. This has the expected 
1/fc^ dependence on large scales, characteristic of a biased tracer of the matter field. We project 
how well upcoming 2i cm observations will be able to disentangle this signal from foreground con- 
tamination. Although foreground cleaning inevitably removes the large-scale modes most impacted 
by primordial non-gaussianity, we find that primordial non-gaussianity can be separated from fore- 
ground contamination for a narrow range of length scales. In principle, futuristic redshifted 21 cm 
surveys may allow constraints competitive with Planck. 



I. INTRODUCTION 

The standard cosmological model makes predictions 
that have been confirmed in a variety of arenas, from the 
cosmic microwave background to large surveys of galax- 
ies. The model contains three mysteries, elements that 
are not part of the Standard Model of particle physics: 
dark matter, dark energy, and inflation. Inflation cur- 
rently plays the role of providing the seeds of structure, 
and the simplest inflationary models predict that the per- 
turbations responsible for this structure were drawn from 
a gaussian distribution. Evidence for primordial non- 
gaussianity then would speak to either a more complex 
inflation model or an alternative in which early acceler- 
ation does not occur. Either would be fascinating and 
probe physics operating at the earliest moments in the 
history of our Universe. 

The Planck Collaboration [l| has, however, recently - 
as we were completing this work - placed stringent con- 
straints on primordial non-gaussianity, by determining 
a robust upper limit to the 3-point function of the cos- 
mic microwave background [3tJ1- These results provide 
support for the simplest models of inflation with gaus- 
sian primordial fluctuations; nature may not provide us 
with this potential handle on the physics of inflation and 
alternatives. Given the enormous significance of a detec- 
tion of primordial non-gaussianity, it is nonetheless worth 
pursuing additional observational constraints. Since the 
cosmic microwave background is relatively unscathed by 
gravitational, non-linear effects, it will be challenging to 
improve on the Planck constraint or to confirm it in- 
dependently. The two-point function of biased tracers 



(galaxies, clusters, etc.) p may, however, be able to pro- 
vide competitive constraints. 

While galaxy surveys have been the premier method 
of studying the large scale structure of the Universe, up- 
coming 21 cm surveys will map the distribution of neu- 
tral hydrogen and therefore provide another handle. Ul- 
timately the large number of Fourier modes potentially 
accessible to 21 cm surveys, may allow even more strin- 
gent constraints than possible with galaxy surveys and 
the cosmic microwave background [y]. In principle, an 
extremely futuristic redshifted 21 cm survey could even 
detect the level of non-gaussianity expected from the sim- 
plest single field models of slow-roll inflation [7|. The first 
aim of redshifted 21 cm surveys is, however, to map out 
the details of the reionization process at z ~ 6 — 20. Mo- 
tivated in part by the ultimate promise of the redshifted 
21 cm line as a probe of primordial non-gaussianity, we 
examine here whether measurements during the reion- 
ization epoch might themselves provide a useful probe. 
During reionization, initial halo formation triggered early 
star formation that produced radiation sufficient to ion- 
ize large bubbles. An important realization from the past 
decade of theoretical work is that regions with large scale 
over-densities are ionized before typical regions (e.g., [8|- 
[ll|V This results because small scale halo formation 
is biased: halos preferentially form in large scale over- 
densities. Given that primordial non-gaussianity affects 
this biasing [5|, it is interesting to study the impact of 
non-gaussianity on reionization. 

There have been several related studies in t he p ast |12| - 
Il5j . Our work has most overlap with Ref. [l^- These 
authors focused on the scale-dependent clustering of the 



ionized regions, while we additionally quantify the impact 
of primordial non-gaussianity on the timing of reioniza- 
tion and the size distribution of ionized regions, and com- 
pare with analytic predictions for the scale-dependent bi- 
asing signature. Furthermore, we consider the impact of 
foreground cleaning in more detail than previous authors. 
This may ultimately provide the largest obstacle for ob- 
taining precise constraints on primordial non-gaussianity 
from redshifted 21 cm observations, and so we estimate 
its impact carefully here. 

The outline of this paper is as follows. We study non- 
gaussianity and reionization here by combining the tech- 
nologies developed by Furlanetto et al. Q to model reion- 
ization and the path integral formalism of Ref. [la] to 
quantify the impact of non-gaussianity. In ^TTl we review 
the Furlanetto model and explain how non-gaussianity 
impacts the collapse fraction and therefore the over- 
density threshold above which a large scale region will be 
considered ionized. Then, we carry out semi-numerical 
simulations in i jIIII and show results for the bubble size 
and ionization history in non-gaussian models. Then, in 
i^IVl we compute the two-point function of the ionization 
field in the simulations and demonstrate that the power 
spectrum rises on large scales. Finally, we conclude in 
§V]by projecting how well upcoming 21 cm surveys will 
be able to measure this feature in the power spectrum 
in the presence of astrophysical foregrounds that pollute 
the large scale spectrum. We conclude in WIl 



II. THE FZH MODEL AND 
NON-GAUSSIANITY 

In effort to understand the impact of primordial non- 
gaussianity on the reionization process, we will extend 
the analytic reionization model of Ref. [8| (hereafter 
'FZH') to the case of non-gaussian initial conditions. In 
order for this work to be self-contained we briefly sum- 
marize the FZH model here, but refer the reader to the 
original paper for a complete treatment. The crux of FZH 
is that galaxies form first in large scale overdense regions 
and that these regions hence reionize before typical parts 
of the Universe. Extended Press-Schechter (EPS) the- 



ory and the excursion set formalism |17l | describe how 
halo formation - and by extension galaxy formation - 
is enhanced in large scale overdense regions, and so we 
can apply EPS to model reionization. This technique 
most faithfully captures large scale variations in the tim- 
ing of reionization; we anticipate that primordial non- 
gaussianity will have its most dramatic impact on pre- 
cisely these large scale variations. Furthermore, upcom- 
ing redshifted 21 cm surveys will measure only the large 
scale features of reionization (e.g. [1^, [I^). For these 
reasons, the FZH model is well-suited for our present 
purposes. 

In order for a region to be reionized, the number of 
photons emitted by sources in the region must at least 
exceed the number of hydrogen atoms contained within 



the region. Since ionized atoms can recombine, it takes 
more than one photon per atom to reionize a patch of 
the Universe: recent simulations suggest that a few pho- 
tons per atom should suffice (e.g. [20|). In the simplest 
variant of FZH that we follow here, one supposes that 
each galaxy, with total host halo mass Afgai, can ionize a 
mass of hydrogen (accounting for some average number 
of recombinations) proportional to its host halo mass. 



Mi, 



CMg 



(1) 



Here C is an ionizing efficiency parameter that depends 
on the fraction of galactic baryons that arc converted into 
stars, the number of ionizing photons that are produced 
per baryon converted into stars, the fraction of ionizing 
photons that escape galactic host halos and make it into 
the surrounding IGM, the average recombination rate in 
the IGM, and other factors. Plausible values for these 
quantities yie ld C '--^ 10, although with substantial uncer- 
tainties (e.g. |18j). 

We further assume that every dark matter halo above 
some minimum (total, i.e., dark matter plus baryons) 
mass, Mmin, hosts a galaxy. Throughout we assume that 
Mmin is set by the mass scale at which the virial tem- 
perature is 10* K, above which gas can dissipate thermal 
energy by emitting atomic lines, condense into the cen- 
ter of the halo, and eventually form stars. It provides a 
plausible lower limit for the host halo mass of a galaxy. 
At high redshifts where Qmiz) ~ 1, the corresponding 
atomic cooling mass scale is j2l| : 



M„,in = 1.2 X lO'^M, 




3/2 ^Q^27N 1/2 



fin 



(2) 



Given these assumptions about the ionizing sources 
and recombinations in the IGM, FZH consider spheres 
of varying radius around every point in the IGM. Each 
gas parcel in the IGM is approximated to be either com- 
pletely ionized or completely neutral. According to Equa- 
tions [1] and m if a sufficiently large fraction of the mat- 
ter contained in a given sphere is collapsed into galaxy- 
hosting halos, the region should be ionized by the sources 
within. In particular, suppose a region has overdensity 
dm when the linear density field is smoothed with a (real 
space) spherical top-hat of co- moving radius _R,„ , enclos- 
ing a Lagrangian mass Mm ~ (pm)47ri?^/3, with (p,n) 
indicating the co-moving cosmic mean matter density. 
Let us further denote the variance of the linear density 
field on this smoothing scale by Smi and the variance 
smoothed on mass scale M^-^ by Sn- The condition for 
the region to be ionized is then: 



C/coU ("Sn I ^m,S'm) > 1- 



(3) 



Here the symbol /coii refers to the collapse fraction in 
the region, i.e., to the fraction of matter in the region 
that is in halos of mass larger than Mmin- The minimum 
overdensity for which this condition is satisfied is given 



the label (5m = ^x- A given point in the IGM is consid- 
ered to be part of an ionized region of size Rm when R^ 
is the largest smoothing scale for which this condition is 
satisfied.^ In the event that there is no smoothing scale 
around a given point for which this criterion is satisfied, 
the point is considered to be completely neutral. 



With these assumptions, FZH treat reionization - in 
the language of excursion set theory - as a 'barrier- 
crossing' problem. In the excursion set formalism one 
considers the behavior of the smoothed field about a 
point, Sm, as a function of decreasing smoothing scale or 
equivalently increasing S„i; each realization of Sm with 
increasing variance, Sm, is said to follow a 'trajectory'. 
The statistical properties of the ionized regions then fol- 
low from considering the probability distribution that 
trajectories cross the barrier condition of Equation [3] at 
various smoothing scales. 



In order to extend this treatment to the case of non- 
gaussian initial conditions, two steps are hence involved. 
First, we need to calculate the conditional collapse frac- 
tion in an overdense region, /coii(5'ri|i5TO, S'm), for non- 
gaussian models. This in turn defines the reionization 
barrier - Sx{Sm) ~ through Equation[3l Second, we need 
to consider non-gaussian modifications to the probabil- 
ity distribution for trajectories to cross this barrier. A 
technical challenge with both of these steps is that non- 
gaussianity induces distinctive correlations between Sm at 
different smoothing scales, which must be incorporated 
into these calculations. This complicates things consid- 
erably in comparison to the case considered by FZH. In 
the usual FZH model, barrier crossing probability distri- 
butions are calculated using a top-hat smoothing filter in 
/c-space and assuming gaussian initial conditions.^ With 
these assumptions, different steps in a trajectory are un- 
correlated: the future evolution of a trajectory with in- 
creasing variance is independent of its past history. This 
is the defining characteristic of a Markov process. For- 
tunately, Maggiore, Riotto and De Simone QQ, |22| - [25| 
developed a path integral formulation of excursion set 
theory which allows one to study departures from the 
Markov case, including the mode- couplings induced by 
primordial non-gaussianity (see also the related works 



FZH describe how this approximately accounts for the possibihty 
that a region is ionized by a neighboring cluster of sources. 
Note that in FZH, and many other applications of excursion set 
theory, formulas for the collapse fraction and other quantities of 
interest are derived assuming a top-hat filter in fc-space, yet ap- 
plied using a top-hat in real space. More specifically, the collapse 
fraction formulas involve the variance of the linear density field 
smoothed on various scales, and these are generally calculated 
using a real space top-hat although the formulas themselves are 
derived using fc-space filters. 



A. Non-gaussian Models 



Throughout the present work, we specify to the special 
case of local models of primordial non-gaussianity. In 
these models, the primordial curvature perturbation (on 
scales smaller than the horizon) is parametrized by: 



<P{X) = <i>G[x) + /nL [4>Ux) - {4>Ux))] 



(4) 



Here (jioix) is a gaussian random field and /nl charac- 
terizes the strength of the non-gaussianity. This form 
produces a bispectrum that is peaked for squeezed trian- 
gles - i.e., for /c-space triangles in which one wavevector 
has much smaller magnitude than the other two. Al- 
though the above form of primordial non-gaussianity is 
only one of many possibilities, it is the most well-studied, 
and is expected for a wide range of different scenarios, 
such as multi-field inflationary models. Recent Planck 
results provide tight constraints on these models, finding 
/j|^°£ai ^ 2.7 ± 5.8 at 68% confidence level Q. We will 
nonetheless consider significantly larger values for /nl 
in order to best illustrate the effects of primordial non- 
gaussianity on reionization. In addition to the local-type 
non-gaussianity considered here, it is also common to 
consider equilateral type non-gaussianity, which peaks for 
triangles with fci « fc2 ~ /ca, as well 'folded' and 'orthog- 
nal' triangle configurations (see e.g. [l[). It may also be 
possible to improve on the Planck coUabroation's current 
68% confidence constraints on these other types of non- 
gaussianity, /n1"" = 42 ± 75, and f^l^° = -25 ± 39 d, 
but we don't consider this explicitly here (although see 



The Reionization Barrier in Non-gaussian 
Models 



We first consider how the reionization barrier is modi- 
fied in models with primordial non-gaussianity. Adshead 
et al. [201 and D'Aloisio et al. [3l| calculate the condi- 
tional collapse fraction in models with primordial non- 
gaussianity using a top-hat filter in k-space and the path 
integral formulation of the excursion set formalism from 
\W{ . For the special case of spherical collapse, Adshead 
et al.'s Equation (46) [30| gives an approximate form for 
the collapse fraction in a region of large-scale ovcrdcn- 
sity Sm- Retaining a few terms, dropped in the large 
scale limit that was taken in Adshead et al. Equation 



(46), wc have: 

/coll(,'-'n|Om: "-"i 




{{€5. 



2{SlSn)) 



Sciz) 



Sm v/27r(5„ - Sraf 



X cxp 






(5) 



Here 5c{z) = 1.686£'(0)/£'(z) is the critical overden- 
sity for spherical collapse scaled to z = and D{z) is 
the linear growth factor. The quantity S'„ is the vari- 
ance of the linear density field, when the density field is 
smoothed on the scale Mmin, Sm is the linear variance 
smoothed on large scale (Mm), and (J^) is the skewness 
at smoothing scale Mmin- Similarly, (S^) is the skew- 
ness smoothed on the large scale Mm, while (S^Sm) and 
(Sm^n) are cross terms that also arise in non-gaussian 
models. Each of these quantities is linearly evolved to 
the present day (z = 0). The first term is the usual re- 
sult for gaussian initial conditions from Press-Schechter 
theory, while the second and third terms are approxi- 
mate modifications for a Universe with non-gaussian ini- 
tial conditions. As detailed in Ref. [30|, the non-gaussian 
corrections - i.e., the second and third terms above - 
are derived assuming (5^(z) ^ Sm and including only 
dominant three-point correlator terms. This expression 
is equivalent to D'Aloisio et al.'s Equation (32) [3l| in 
the limit S'^{z) ^ Sm, Sc{z) ^ 6m, and ignoring their 
last term (i.e., the last line of their Equation (32), pro- 
portional to their 'C), which is negligible for the case 
considered here. 

For positive /nl, the conditional collapse fraction of 
Equation [5] is enhanced compared to the case of gaus- 
sian initial conditions. This results because positive /nl 
enhances the high density tail of the probability distri- 
bution function and thereby boosts the chance of be- 
ing above the collapse threshold, Sc{z). In the large 
scale limit, where 5'„ ^ Sm, the correlators (S^) and 
{SmSn) are small compared to (Sf^) and (Sf^Sm), and 
so the corresponding terms in Equation [5] are only sig- 
nificant close to the minimum mass scale, Mmin (Sn)- 
Furthermore, note that in the very large scale limit, 
(S'^Sm) / Sm oc \J\I Sm ~^ oo, i.c., this term diverges to- 
wards large scales."^ This divergence is not a concern 
since the density contrast, 6m, is itself tending towards 



^ We have assumed the Newtonian form of the gravitational poten- 
tial and so our calculation breaks down on scales near the hori- 
zon, where relativistic effects are important. One might worry 
that this divergence is an artifact of assuming the Newtonian 
form for the gravitational potential. Wands & Slosar [32l |. how- 
ever, carry out a relativistic analysis and find that the halo bias 



zero on large scales. The (S^Sm) correlator in the expres- 
sion here describes the mode-coupling between large and 
small scales and it is this term that ultimately leads to 
the scale-dependent halo clustering signature (e.g., J3G|'). 
Presently, we are interested in the impact of this term on 
the reionization barrier. Since it causes the conditional 
collapse fraction to blow-up on large smoothing scales 
(for positive /nl), it leads to a down-turn in the reioniza- 
tion barrier at small Sm- In the case of negative /nl, the 
sign of this effect is reversed and the reionization barrier 
turns-up at small Sm- We caution that approximations 
made in Equation [5] may influence the shape of the bar- 
rier on the largest scales here, and so we caution against 
taking this too literally. The behavior of the barrier on 
these scales does not, however, impact our results. 
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FIG. 1: FZH reionization barrier for various /nl models. The 
lines show the critical overdensity for a region to be ionized 
as a function of the linear variance when the density field 
is smoothed on the scale of the region (with both quantities 
linearly extrapolated to 2: = 0). Each curve is for z = 9, 
and assumes that galaxies form in halos above the atomic 
cooling mass, Mcooi, with an ionizing efficiency parameter of 
C^ — 12. For /nl = 0, the volume- averaged ionization fraction 
is (xi) — 0.48 at this redshift. The collapse fraction in a 
region of large-scale overdensity increases with increasing /nl 
and so the critical overdensity for a region to be ionized, 5x 
decreases with increasing /nl. The black dot-dashed line is a 
/nl = model with C adjusted to match the volume-averaged 
ionization fraction in the /nl = 300 model (see text). 



(and hence the conditional collapse fraction considered here) in- 
deed diverges on large scales. 



In order to quantify the impact of /nl on the shape of 
the reionization barrier, we invert Equations [3] and [5] to 
find 6x as a function of Sm for several /nl models. Here 
we generally consider 2 = 9, C = 12, and set Mmjn to 
the atomic cooling mass corresponding to a virial tem- 
perature of Tvir = W^K, as specified by Equation [2j For 
/nl = these parameters give (xi) ~ 0.48, and so we 
are considering roughly the 'mid-point' of reionization, 
where half of the volume of the Universe is ionized. The 
results of these calculations are shown in Figure [1] for 
models with /nl = (-300,-100,0,100,300). The bar- 
rier decreases with increasing /nl: positive /nl boosts 
the conditional collapse fraction in overdense regions, and 
hence a lower critical overdensity Sx is required for a re- 
gion to be ionized in a non-gaussian model. The barrier 
turns down (up) significantly on large scales for positive 
(negative) /nl because the collapse fraction blows up to- 
wards large scales, as mentioned earlier. Even for values 
as large as |/nl| ~ 100 - now strongly disfavored by the 
Planck constraints [l[ - this down-turn (up-turn) occurs 
on rather large scales, where the barrier-crossing proba- 
bility should be very small. In general, /nl leads to only 
small changes in the barrier height and shape. For in- 
stance, Sx is ^ 8% smaller at _R = 3 Mpc/ft, {Sm ~ 3) 
and -- 15% smaller at i? = 10 Mpc/h (5,„ « 0.5) in the 
/nl = 100 model compared to the gaussian, /nl = 0, 
case. The scale i? = 3 Mpc/h mentioned here corre- 
sponds to the peak in the analytic bubble size distribu- 
tion for the /nl — model, computed as in [8|, while 
-- 99% of bubbles in this model have R < 10 Mpc//i, 
the second scale considered here. These numbers hence 
give some indication of which scales typically cross the 
reionization barrier in these models. 

In addition, the impact of /nl on the reionization bar- 
rier at a given redshift is largely degenerate with the 
effect of varying (. This is illustrated by the black 
dot-dashed line in Figure [TJ which shows an /nl = 
model with ( enhanced (from ^ = 12 to C = 14.3) 
to match the volume-averaged ionization fraction in the 
/nl ~ 300 model at this redshift. The shapes of the bar- 
riers are somewhat different; the barrier in the /nl model 
has the large scale down-turn and rises more steeply on 
small smoothing scales. The non-gaussian model should 
have slightly more large bubbles and slightly fewer small 
bubbles than a gaussian model with the same volume- 
averaged ionization fraction (see ^IIII for further details) . 
However, from a similar bubble-size distribution calcu- 
lation to the one mentioned above, we expect ^ 99% of 
random walks in the /nl = 0, C = 14.3 model (that is 
largely degenerate with the /nl — 300 model) to cross 
the barrier at Sm J^ 0.4. This suggests that the most 
prominent differences between the two barriers, on large 
smoothing scales, will have little impact on the resulting 
bubble-size distributions since few random walks cross 
the barrier on such large scales, at least near the mid- 
dle of reionization. Hence it seems that varying ^ should 
largely compensate for /nl induced changes in the reion- 
ization barrier. In addition, we should keep in mind that 



the values of /nl considered here are already strongly dis- 
favored by existing data, and so the degeneracy is even 
more important than in this illustrative case. 
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FIG. 2: Volume-averaged ionization fraction vs. redshift for 
various /nl models. The curves show the impact of /nl on the 
reionization history of the Universe in models with Mmin = 
Mcooi and ^ = 12. The reionization process is accelerated in 
a Universe with positive /nl, while it is delayed in a Universe 
with negative /nl- 



A similar calculation determines the volume-averaged 
ionization fraction as a function of redshift in the FZH 
model. Specifically, this is given by: 



Cerfc 



C/coll(<S'n) 

Sc{z) 



+ c- 



{sf.) 



3V27rS', 



3/2 






-1 



exp 



25„ 
(6) 



This equation holds for redshifts above which the ion- 
ized fraction, (xi), becomes unity. Here fcon{Sn) is the 
(non-conditional) collapse fraction for halos above the 
minimum mass at any given redshift. This is specified 
by Equation [5] in the limit that Mm — > 00. Sm — > 0. 
Since the collapse fraction is enhanced in models with 
positive /nl, so is the volume-averaged ionization frac- 
tion. This is quantified in Figure [U which shows the 
volume-averaged ionization fraction as a function of red- 
shift in the /nl models of Figure [H As before, we fix 
C = 12 in each model and the minimum host halo mass 
at the atomic cooling mass. Reionization starts and fin- 
ishes earlier (later) in models with positive (negative) 
/nl compared to models with gaussian initial conditions. 
However the effects are small: for instance, the ioniza- 



tion fraction near [xi 
hi. = 100. 



0.5 is boosted by < 10% for 
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FIG. 3: Contours of constant ionized fraction in the (" — /nl 
plane. The blue lines are lines of constant ionized fraction, 
with {xi{z = 9)) = 0.3,0.5,0.7,0.9 and 0.99 respectively, il- 
lustrating the degeneracy between /nl and C,. 



The impact of /nl on the volume-averaged ioniza- 
tion fraction is, however, degenerate with uncertainties in 
the ionizing efficiency and the minimum host halo mass, 
which are unknown and in any case provide only a rough 
model for {xi){z). We illustrate this degeneracy in Fig- 
ure|3l which shows contours of constant {xi{z = 9)). The 
plot spans a very large range in /nl, including models 
that are already strongly disfavored by existing data for 
illustration. The steepness of the contours indicates that 
relatively small variations in C, can compensate for /nl- 
induced changes, even over the large range in /nl shown. 
In Figure [T] we found that the reionization barrier in an 
/nl model nearly matches that in a gaussian model, once 
we adjust C, to fix (2;^) across models. This suggests that 
the bubble size distribution will mostly share this degen- 
eracy: in other words, the contours of fixed {xi) in Figure 
[3] should resemble contours of fixed bubble-size distribu- 
tion as well. We will explore more distinctive and unique 
imprints of primordial non-gaussianity subsequently. 



III. SEMI-NUMERIC SIMULATIONS AND 
NON-GAUSSIANITY 



tion [3] at various smoothing scales. Here we will handle 
this numerically, using the so-called 'semi- numeric' sim- 
ulation technique developed in Q ■ This simulation tech- 
nique is essentially a (three dimensional) Monte-Carlo 
implementation of the FZH model. The Monte-Carlo im- 
plementation has the advantage that it partly captures 
asphericity in the shapes of the ionized regions, (which 
we will sometimes refer to as 'ionized bubbles'). Further- 
more, it provides mock reionization data cubes that are 
convenient for measuring the statistical properties of the 
epoch of reionization, and for visualizations. In compar- 
ison to radiative transfer simulations of reionization, the 
semi-numeric [9| technique has the advantage that it is 
extremely fast, while still capturing the large scale fea- 
tures of the reionization process fairly accurately [3, [Sj] . 
In order to use this technique to study reionization in 
an /nl model, we need to first generate a non-gaussian 
realization of the linear density field in the model of in- 
terest. We do this in the usual manner, briefly described 
here for completeness. Specifically, we start by generat- 
ing a gaussian random realization of the gaussian field 
(J)q{x). In generating this realization, we assume that 
(1)q{x) has a scale invariant power spectrum of the form 
A2(yfc) = k^P^{k)/{2TT^) = 8.71 X 10-1°.-* From (t>G{x), 
we generate the primordial curvature perturbation as- 
suming the local model described by Equation 21 The 
density field then follows from the potential perturba- 
tion by Poisson's Equation, which may be written (and 
applied) in Fourier space as: 
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Finally, we Fourier-transform the resulting density field 
into real space. In each model considered, the same set 
of random numbers are used to generate the underly- 
ing gaussian random potential field, (J)g{x). This lessens 
the impact of sample variance when comparing differ- 
ent models, and isolates the impact of primordial non- 
gaussianity. 

We generate models with /nl = 

-100,-50,0,50,100,300 in two different simulation 
volumes. The first simulation volume has a co-moving 
side length of Lbox = 150 Mpc//i, while the second simu- 
lation has Lbox = 2 Gpc//i. In each case, the density and 
ionization fields are tabulated on a 512'^ Cartesian grid. 
The smaller simulation box captures small reionization 
bubbles, while the larger volume runs are essential 
for examining the large scale clustering of the ion- 
ized regions, especially the distinctive scale-dependent 
signatures induced by primordial non-gaussianity. 

In order to construct the ionization field from a real- 
ization of the linear density field, we use the procedure 
of FZH and Q- We smooth the linear density field on 



In order to calculate the statistics of the ionized re- 
gions, such as their size distribution, we need to consider 
the probability that trajectories cross the barrier of Equa- 



* For reference, this model has eg : 
recent constraints, e.g. |34|| . 



0.86, broadly consistent with 



a range of scales, starting from large scales and gradu- 
ally stepping down to the size of the simulation pixels. 
A pixel is marked as ionized if it crosses the barrier of 
Equations [3] and [5] on some smoothing scale, while pix- 
els that fail to cross the barrier on any smoothing scale 
are marked neutral. Since we have proper non-gaussian 
realizations of the density field, the enhanced probabil- 
ity of crossing the rcionization barrier in /nl raodels is 
naturally accounted for in this step. 
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FIG. 4: Slices through semi-numeric simulations of reioniza- 
tion models with primordial non-gaussianity. In each model, 
(" = 12, z = 8. The slices are 150 co-moving Mpc/h on a side, 
and are 0.29 Mpc/h thick. Top Left: /nl = -100, Top Right: 
/nl = 0, Bottom Left: /nl = 100, Bottom Right: /nl = 300. 
The black regions show ionized bubbles, while the red regions 
are neutral gas. As /nl increases, reionization is more pro- 
gressed (for fixed ionizing efficiency, (, and redshift, z), and 
the bubbles are larger. 



We show example slices through 150 Mpc/h semi- 
numeric reionization simulations with ^ = 12 and z = 8 
in Figure ID Since each simulation is generated with the 
same underlying gaussian random part of the potential 
field, (J)g{x), we can compare the slices directly, region- 
by-region. As expected, reionization has progressed fur- 
ther, and the ionized regions are larger, in the models 
with primordial non-gaussianity. In the /nl ~ ±100 
model, however, the differences with the case of gaus- 
sian initial conditions are somewhat subtle. Note that 
we are comparing the different models at fixed C, and z, 
and so the models have varying volume-averaged ioniza- 
tion fractions, which increase with /nl- Specifically, the 
model with /nl = has {xi) = 0.45 at z = 8, while 
{x,) = (0.40,0.50,0.63) for /nl = (-100,100,300) at 
z ~ 8. A caveat here is that these values are smaller 
than expected from the analytic curves in Figure [2l and 
the difference with the expected value decreases with in- 
creasing |/nl|- The smaller values result for two reasons: 
first, there are some ionized bubbles that are smaller 
than the size of our simulation grid and hence not cap- 
tured. This effect is more important at small |/nl| since 
the bubbles are smaller in these models. Second, our 
collapse-fraction expressions are derived using a top-hat 
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FIG. 5: Bubble-size distribution for various /nl models. The 
curves show the probability distribution of the sizes of the 
ionized regions for xth = 0.9 (see text). The models are at 
z = 8 and most assume (" = 12, Mmi-n ~ Mcooi, for various 
values of /nl as labeled. The black dot-dashed model is for 
an /nl = model with (" adjusted to match the mean ionized 
fraction {{xi)) in the /nl = 300 model. This illustrates that 
the variations in bubble size with /nl at fixed {xi) are much 
smaller than at fixed ^. 



filter in /c-space but applied with a real-space filter: as 
discussed in the Appendix of [9| , this leads to departures 
from the expected global ionized fraction. Since the de- 
parture decreases with increasing /nl, our semi- numeric 
simulations likely overestimate the impact of /nl on the 
ionized fractions and bubble-sizes, but they still serve to 
illustrate the main effects of non-gaussianity on rcioniza- 
tion. 

In order to quantify the visual impressions of Figure 
m we calculate the probability distribution of the sizes 
of the ionized regions for each /nl model. This depends 
somewhat on one's definition for the size of the com- 
plex, aspherical ionized regions in the simulation. Here 
we define the size of the simulated ionized regions as in 
[01 . Briefly, we spherically average the ionization field 
on various scales R, starting from large scales, stepping 
downward in size until wc eventually get to the size of our 
simulation pixels. At each smoothing scale, we compare 
the spherically averaged ionization field to a threshold 
ionization value, xth- A simulation pixel is considered 
'ionized' and belonging to a bubble of radius R, when R 
is the largest smoothing radius at which the smoothed 
ionization field crosses the threshold. Pixels that do not 
cross the threshold on any smoothing scale are consid- 
ered neutral. This procedure is of course similar to the 



way in which the ionization field is constructed in the 
first place. 

The resulting probability distribution function (PDF), 
for xth ~ 0.9, is shown for our various /nl models with 
C = 12, z = 8 in Figure[5l Note that each PDF is normal- 
ized to 1 rather than to the ionized fraction: the PDFs 
show the fraction of bubbles that have a radius between 
R and R + dR. To mention one quantitative description 
of these results, we examine by how much the character- 
istic size of ionized bubbles - which we identify with the 
peak in the bubble size PDFs of Figure [5] - varies with 
/nl for the present values oi C,z. Compared to a model 
with gaussian initial conditions, the characteristic bubble 
size increases by a factor of 1.4 for /nl = 100, a factor 
of 3.7 for /nl = 300, while it decreases by a factor of 1.5 
for /nl = —100. As anticipated in the previous section, 
the bubble-sized distribution is degenerate with (. The 
black dot-dashed line shows that an /nl — model with 
^ adjusted upward to match the ionized fraction in the 
/nl = 300 model. The resulting bubble-size distribu- 
tions are quite similar, although the bubbles are a little 
bit larger in the /nl model. This is expected from Figure 
[I] the main difference between the /nl barrier and the 
gaussian barrier at fixed (xi) is that the /nl barrier is a 
little lower on larger scales (smaller Sm ) allowing slightly 
larger ionized regions to form. 

One might wonder whether random walks can start 
to cross the reionization barrier at very low Sm in non- 
gaussian models, where the barrier has this distinctive 
downturn. In principle, this might lead to a bi-modality 
in the bubble size distribution: perhaps the downturn 
allows some trajectories to cross the barrier on very 
large scales that would be prohibited from crossing oth- 
erwise. On smaller smoothing scales, the trajectory- 
crossing probability might shrink as the smoothing scale 
becomes smaller than the down-turn scale and the bar- 
rier increases, until trajectories catch up again close to 
the scale of the usual peak in the bubble size distribu- 
tion. This might imprint a distinctive large-scale bump 
in the bubble size distribution. It is possible that this 
is even the origin of the kink in the high-i? tail of the 
/nl = 300 model in Figure [SJ This effect might be more 
realizable at the end of reionization, when random walks 
start crossing the reionization barrier on progressively 
larger scales. However, this effect is unlikely important 
for the much smaller values of /nl presently allowed, and 
so we don't investigate it further here. It also may be an 
artifact of approximations made in deriving Equation [SI 
as mentioned earlier. 



IV. SCALE-DEPENDENT BIAS 

The previous results describe the general impact of 
primordial non-gaussianity on reionization. However, 
the most promising approach for obtaining observational 
constraints on non-gaussianity from reionization studies 
is to measure the scale dependent clustering of the ion- 



ized regions. This is directly analogous to the case of the 
clustering of dark matter halos. In the case of halo clus- 
tering, Dalai et al. Q showed that local non-gaussian 
models give rise to a scale-dependent clustering signa- 
ture. If a similar signature arises in the clustering of 
ionized regions, this may provide a distinctive indicator 
and allow constraints on /nl from the epoch of reioniza- 
tion, despite uncertainties in the properties of the ioniz- 
ing sources, and the surrounding IGM. 

Indeed, we would expect the ionized regions to have 
a similar scale-dependent clustering signature to that of 
the dark matter halos. In terms of excursion set model- 
ing, the reionization case is different than halo clustering 
only in that the reionization barrier (Equation |3]) has a 
different shape than the halo collapse barrier. Otherwise 
the physics underlying the scale dependent clustering is 
similar. In particular, a positive value of /nl implies 
that the variance of the density field is enhanced in large 
scale regions with above average potential perturbation 
(which are overdense on large scales as quantified by the 
Poisson Equation, Equation [T]) The enhanced variance 
boosts the collapse fraction in such regions, and increases 
the tendency for these regions to be ionized before typ- 
ical regions. In the context of the FZH model, the ex- 
tra variance in large scale overdense regions boosts the 
rate at which trajectories cross the barrier of Equation [3] 
just as it increases the rate of crossing the spherical col- 
lapse barrier in the case of halo clustering. In both cases 
the form of Equation [7] then implies a distinctive scale- 
dependent clustering term, ex l/{k'^T{k)) Q. Moreover, 
we [301 (see also [3l[) showed that this general form is 
expected for an arbitrary barrier crossing problem in the 
presence of primordial non-gaussianity with a non-zero 
squeezed limit. One caveat in the reionization case is 
that we expect this to apply only on scales much larger 
than the size of the ionized regions. On smaller scales, 
the ionization field will decorrelate from the underlying 
density field since small scale regions are either highly 
ionized or completely neutral irrespective of the precise 
value of the density field. 

Our aim here is to quantify the scale-dependent clus- 
tering of the ionized regions. In order to capture the large 
scale Fourier modes where the scale-dependent cluster- 
ing signature should dominate, we use 2 Gpc/h semi- 
numeric reionization simulations. As before, we con- 
struct simulations for several values of /nl; here we 
simulate /nl = 0, ±50, ±100. We further calculate the 
power spectrum of the ionization field in each model, 
Px,x{k).^ We also calculate the ionization-density cross 
power spectrum, Px,s^{k) and the density power spec- 
trum, Ps ,5 (k). We then estimate the bias of the ion- 
ization field from the ratio of the ionization-density cross 



^ We consider the power spectrum of Xi rather than Sx = {xi 
{xi))/{xi), i.e., we do not normaUze by (xi). 



power spectrum and the density power spectrum, 



(8) 



As discussed previously, we expect primordial non- 
gaussianity to slightly lower the collapse barrier of Equa- 
tion[3](for positive /nlj as illustrated in Figure[T]), and for 
the enhanced high density tail in these models to increase 
the probability of crossing the reionization barrier. The 
latter effect, and the coupling between large and small 
scale modes, is responsible for the scale dependent clus- 
tering enhancements. However, each of these two effects 
modifies the bias of the ionized regions, hx{k). Since /nl 
is known to be small, we can consider each effect as a 
small correction in a Taylor expansion around /nl = 0, 
and calculate the two effects separately. The change in 
bias from the reduced barrier height in a positive /nl 
model is mainly due to the fact that the ionization frac- 
tion is larger (than gaussian). For small changes in /nl, 
and on scales much larger than the ionized regions, this 
leads to a roughly scale- independent change in bj.{k). 
We find that this scale-independent term is fairly well 
matched by considering the change in bias in an /nl = 
model after adjusting C, to match the enhanced (a;,;) in 
the corresponding /nl model (see Figure [6] for an illus- 
tration). 

The second effect - increased probability of barrier 
crossing in a non-gaussian model - produces a distinctive 
scale-dependent clustering enhancement of the [5| form: 
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(9) 



Here /S,hx{k) denotes the change in bias with wavenum- 
ber k owing to /nl, and the above equation describes 
only the scale dependent contribution. The quantity 6^ 
is the bias of the ionization field in a gaussian model, and 
T{k) is the transfer function. Here 5b is a proportionality 
constant related to the height of the reionization barrier. 
Note that &^ and 5b depend on both redshift/ionization 
fraction and the reionization model. The growth fac- 
tor here, D{z), is the linear growth factor normalized to 
1/(1 -|- z) during the matter-dominated era. 

Our numerical results mostly perform this perturbative 
analysis, as shown in Figure [6] for z = 8 and ( = 12.^ 
We adjust the single parameter 5b in Equation [S] above 
to match the simulation results, and calculate the scale- 
independent correction to the bias as described above. 
The result is shown by the solid curves in Figure HI The 
fit is a fairly good overall match to the simulation results, 
although some differences are evident. Most importantly. 



The ionized fraction is lower by a factor of ~ 2 than in our 150 
Mpc//i volume because some of the small bubbles are unresolved 
in the large volume simulated here, which focuses on capturing 
the large scale bias. The limited resolution of our calculations 
should impact the bias numbers a little at a given {xi), but should 
not impact the main trends. 
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FIG. 6: Scale-dependent bias of ionization field. The dashed 
lines show measurements of bx{k) (Equation [8]) from semi- 
numeric simulations for z = 8 and fixed C, — 12. The solid 
lines show (by eye) fits to the simulation results. The fits 
consist of two terms. The first term is a scale dependent 
enhancement in the form of Equation [9l with the proportion- 
ality constant, Sb, adjusted to (roughly) match the /nl ~ 50 
measurement. The second term is predicted by varying (" in 
a gaussian model to match the {xi) value in the correspond- 
ing /nl model. As an illustrative example, the dot-dashed 
line shows 6a;(A;) for an /nl ~ model with the same {xi) 
as for /nl = 50. The difference between the dashed and 
dot-dashed lines hence refiects the roughly scale-independent 
(on large scales) correction to the bias for /nl = 50, while 
the enhancement at fc < 0.02h Mpc~^ is dominated by the 
scale-dependent contribution. 



one can see the scale-dependent enhancement of the Dalai 
et al. type on large scales, and that this term is linear 
in /nl, as expected. (Note that the 'fit' curves share 
noisy features with the simulation results, because the 
scale independent correction is calculated directly from 
the simulated power spectra.) 

The fits, however, seem to slightly underproduce the 
large scale enhancement for /nl = 100, and slightly over- 
produce the effect at /nl = —100. The departure from 
the fits at large /nl likely results because the change in 
the scale-independent bias in these models is becoming 
significant and is not perfectly captured by matching to 
a fixed {xi). After all, the reionization barrier and bub- 
ble size distribution still differ slightly after matching to 
a fixed (xi) (see Figure [T]). Along these lines, we find a 
better fit if, in calculating the scale-independent correc- 
tion, we boost C a bit beyond that required to match the 
enhanced (xi) in the non-Gaussian model (although this 
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improvement is not shown in the figure here). Quanti- 
tatively, we find {b^ — 1)Sb ~ 0.49 for the redshift and 
ionization fraction considered here. This coefficient also 
depends on the particular reionization model considered 
here, e.g., the assumption that the ionizing efficiency in 
Equation [T] is itself independent of mass. Note that we 
considered the scale-dependent bias of Xi here, if we had 
instead considered the field S^ = {xi — {xi))/{xi) the scale 
dependent bias coefficient would be ss 2. 

The relatively good match emboldens us to use the 
simple fitting formula of Equation |9] in f|V]to project how 
well future 21 cm surveys will be able to extract /nl- The 
imperfect fit does, however, suggest that more theoretical 
work will be needed to extract the exact value of /nl from 
the data.^ 

Note that the scale dependent enhancement dominates 
over the variation from the enhanced (x) only on very 
large scales (small fc), k < 0.02h Mpc~^. Unfortunately, 
these scales will likely be challenging to observe (see JfV)) . 
In practice, the efficiency of the ionizing sources and their 
host halo masses will not be known a priori and so we will 
likely need to perform a joint fit for the scale indepen- 
dent and scale dependent bias contributions. If {xi){z) 
can be determined from independent observations, this 
would be helpful in separating out the scale dependent 
enhancement. 

We also investigated how the results depend on the 
particular stage of the reionization process, by consider- 
ing an additional model in which ( is enhanced beyond 
our fiducial values of C = 12. In the alternate case con- 
sidered, (xi) = 0.51 at z = 8 (for /nl ~ 0) and so it de- 
scribes the impact of non-gaussianity near reionization's 
midpoint. In this case, we find a similar fit works, espe- 
cially if we allow a slight boost to the scale-independent 
term, beyond that required to match to a fixed (xi). At 
this stage of reionization, we find {b^ — 1)6b ~ 0.35. The 
proportionality constant in the fit of Equation [9] hence 
varies with C,, but this is to be expected since the propor- 
tionality constant should depend on the barrier height, 
which itself depends on (^. In principle, one can also con- 
sider later stages in reionization, but we expect this to be 
less useful. First, the signal drops off towards the end of 
reionization (e.g. [l9[). Second, note that the ionization 
and density fields decorrelate on scales less than the size 
of the ionized regions. As the ionized bubbles grow, only 
larger and larger scales are useful for the scale-dependent 
clustering signature. 

It is also possible to compute the scale-dependent bias 
coefficient, (&^ — 1)(5b, analytically. The analytic calcula- 
tions allow one to quickly explore many /nl models for a 
wide range of ionized-fractions, (xi), and rcdshifts. Ref. 



^ As we were finishing this work, Ref. [3511 was posted, making 
essentially the same point; quantitatively their results may differ 
a bit more from Equation |9] than do ours, but the qualitative 
conclusion of a scale-dependent bias has now been verified by 
several different groups using a range of techniques. 



[30| has derived expressions for the scale-independent 
bias and scale-dependent bias coefficient that arc applica- 
ble to any collapse barrier using a path integral approach 
developed by |36[438| . Following [30|, we define Va{Sm) 
as an expansion in the barrier: 
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The scale-independent part of the bias can then be writ- 
ten as 
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which reduces to the standard expression from e.g. (39[ 
in the limit that the barrier is linear. Since 5x{Sm) and 
Sm are linearly extrapolated to z = 0, the growth factor 
enters here in calculating the bias at redshift z. For an 
/nl cosmology the scale-dependent bias can be written 
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(14) 

In the limit that barrier is fiat, i.e. SxiSm) = ^c, then 
'Poi'Sm) = and the above expression reduces to the 
standard prediction c{Sm)j= Sc{b^ — l)D{z)/D{0), which 
has been derived in e.g. Q. 

To compare these analytic expressions for the bias with 
the results of the semi-numeric simulations presented 
above, we must calculate the volume- weighted average 
of b^{Sm) and c{Sm) over all ionized bubbles. This aver- 
aging is accomplished by integrating over the bubble size 
distribution: 



dmb^{S{m))^V{m) (15) 

dm 

dn 
drac{S{m))— — V{m) (16) 

dm 



where 



{X^) = 



fin 

dm-— Vim) (17) 

dm, 



11 



and V{m) = m/p is the volume of a region of mass m. 
The mass function dn/dm can be computed analytically 
from the excursion set formalism: 



dn 1 p 
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where J^ is the so-called unconditioned crossing rate dis- 
cussed in [301 • 

The results of our calculation of the volume aver- 
aged scale-independent bias and volume average scale- 
dependent bias coefficient are shown in Figures [7] and |H1 
In these figures we have multiplied B^ by {xi{z)) to obtain 
the bias of Xi rather than the field 6x = {xi — {xi))/{xi); 
this allows for direct comparison to the semi-numerical 
simulation results presented above. Figure [7] shows the 
result of calculating the bias for three different collapse 
barriers: a linear fit to the gaussian barrier (black solid 
curve) , the true gaussian barrier (red dashed curve) , and 
the barrier with the effects of non-gaussianity included 
at the level of /nl = 100 (blue dotted curve). The first 
of these barriers corresponds to the calculation of the 
bias presented in @. We see that the effects of non- 
gaussianity at the level of Jnl = 100 on the mean scale- 
independent bias are comparable in magnitude to the 
effects of approximating the barrier as a linear function 
of the variance. 

Figure [5] shows three different approaches for calculat- 
ing the scale-dependent bias coefficient. The black (solid) 
and red (dashed) curves compute this coefficient using 
only the first term of Equation [T3] (i.e. the standard re- 
sult) for a linear barrier and the true barrier, respectively. 
The blue (dotted) curve, on the other hand, includes all 
of the terms in Equation [T31 Including only the first 
term, the linear and true barrier calculations give very 
similar results. Including all of the terms, however, gives 
a significantly larger scale-dependent bias; this difference 
grows with decreasing redshift as reionization proceeds. 
We find that the first term approximation gives better 
agreement with the simulation results. In [30[, we found 
related differences between the first term approximation 
and including all of terms for the case of an ellipsoidal 
collapse barrier (see Figure 3 of [SO]). In that paper, it 
was speculated that approximations used in deriving Eq. 
[13] - in particular, approximating 3-point functions by 
their end-point values - may lead to an artificial rise in 
the scale-dependent bias coefficient when all of the terms 
in Equation 1131 are included. A similar effect may make 
this coefficient spuriously large here. 

The analytic calculations roughly agree with our nu- 
merical results, provided the first-term approximation is 
more robust and so we compare with it and adopt this ap- 
proximation in what follows. The scale-independent bias 
is roughly 1.5 and the scale dependent bias coefficient - 
denoted by c{xi) in the analytic calculation of Equation 
[T3l " is roughly (6^ — 1)5^ sa 0.5 over the range of z that 
we consider. This is close to the simulated values at sim- 
ilar ionized fractions (see Figure [5]). There are several 
reasons why the analytic calculation might not agree ex- 



0.5 0.4 



<x,> 

0.3 



A 
X 
V 1.5 



II 1 1 


-^ 


1^^^ f^j_ = 0, true barrier 


1,1, 



FIG. 7: The volume averaged scale-independent bias. 



actly with the simulation results. For one, the analytic 
calculation assumes that the ionized regions are spheri- 
cal, while it is clear from Figure[3|that this is not the case. 
Second, the analytic results are based on the first-term 
approximation to Equation 1131 In order to robustly cal- 
culate the corrections here, it may be necessary to move 
beyond approximating the 3-point functions at their end- 
point values, as mentioned above. Next, the simulations 
here do not resolve the smallest bubbles, and do not pro- 
duce precisely the expected (xi), as discussed previously. 
Asides for these caveats, we find the agreement between 
the analytic and numerical calculations encouraging and 
will therefore use the analytic calculations to estimate 
the prospects for constraining /nl with future surveys 

Our semi-numeric results regarding scale-dependent 
biasing also agree broadly with Joudaki et al. p^ . 
who first identified the scale-dependent clustering en- 
hancement in the ionization field. These authors ap- 
proximated, however, the reionization barrier as fixed, 
while we further included the impact of primordial non- 
gaussianity on the reionization barrier itself. As a result, 
they did not include the scale-independent enhancement 
discussed here. However, we have verified that this is 
not a big obstacle, and that the scale-dependent signa- 
ture can still be easily discerned. We expect this to be 
even more the case for smaller values of /nl than consid- 
ered here; given the tight constraints from Planck data, 
smaller values of /nl ^ ±10 span the currently interest- 
ing regime. 
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the EoR (e.g. [44| - |46( V° With these approximations, the 
brightness temperature contrast is: 



0.5 

ar-r- 



First term, linear barrier 
First term, true barrier 
All terms, true barrier 



FIG. 8: The volume averaged coefficient of the scale- 
dependent bias. 



IS IT MEASURABLE? 



Now that we have characterized the impact of primor- 
dial non-gaussianity on reionization and quantified its 
signature in the scale-dependent clustering of ionized re- 
gions, we turn to consider whether these signatures may 
be observable in upcoming reionization surveys. The 
most promising approach is to use the scale-dependent 
bias, since the other effects are likely degenerate with un- 
certainties in the properties of the ionizing sources and 
the surrounding IGM. Three main types of observations 
have been discussed in the literature that can poten- 
tially probe or constrain the scale-dependent clustering of 
the ionized regions: narrow-band surveys for Ly-a emit- 
ting galaxies during reionization (e.g. [40,|41|), measure- 
ments of the small-scale CMB anisotropics induced by the 
patchy kinetic Sunyaev-Zel'dovich effect (e.g. [43. l43|). 
and measurements of the redshifted 21 cm line from the 
epoch of reionization (EoR) (e.g. [l^l)- The latter mea- 
surement is the most direct probe of spatial fluctuations 
in the ionized fraction during reionization, and so we fo- 
cus on these measurements here. 

The main quantity of interest for the redshifted 21 cm 
measurements is the 21 cm brightness temperature con- 
trast between a neutral hydrogen cloud at redshift z and 
the CMB. We work in the limit that the spin temper- 
ature of the 21 cm line is much larger than the CMB 
temperature throughout all space, and further, we ignore 
rcdshift-space distortions from peculiar velocities. These 
are expected to be good approximations during most of 



Ti\ = T^XYaiX + (5p). 



(19) 



The constant Tq is Tq = 28 [(1 + z)l\G\l'^ niK for our 
cosmological parameters [T^. Here xhi is the neutral 
fraction, and 8p is the density contrast. Each of T21, 
xhi, and bp vary spatially and with redshift, but we have 
suppressed this dependence in our notation here. 

The power spectrum of 21 cm brightness temperature 
fluctuations is related to the spatial fluctuations in the 
ionization field we considered earlier. If we expand to first 
order in fluctuations in the density contrast and to first 
order in neutral fraction fluctuations (neglecting redshift 
space distortion and spin temperature fluctuations), we 
expect: 



^21,21 ('^■) 



T2[62_2(1-(x,))5, + (1-(x,))2 



xAL,x(fc). 



(20) 



Here hx is defined as in Equation |SJ it is the linear bias 
factor of the ionization field, rather than the bias of the 
fluctuations in the ionization field that is sometimes con- 
sidered. The bias factor wc consider here is equal to the 
bias factor of the ionization fluctuations, multiplied by a 
factor of (xi). In general, working to first order in Sx can 
be problematic since the ionized regions are expected to 
be large during most of the EoR. The large spatial fluc- 
tuations in the ionization fraction imply that additional 
terms, dropped in Equation l20[ can be important even on 
large spatial scales [43] . Since we are interested here only 
in the very large scale clustering signature, that applies 
on scales much larger than the size of the ionized regions, 
the approxiination of Equation ^U\ should nonetheless be 
adequate. In this approximation, Eauation l20l relates the 
bias factor of the 21 cm fluctuations to the bias factor 
measured in the previous section. 

A challenge for measuring the scale-dependent cluster- 
ing signature with future 21 cm measurements is that 
primordial non-gaussianity significantly impacts cluster- 
ing only on rather large scales. On these large scales 
foreground contamination in the redshifted 21 cm data 
may be prohibitive. At the frequencies of interest, fore- 
ground contamination from galactic emission and extra- 
galactic point sources is expected to have a mean bright- 
ness temperature that is roughly four orders of magni- 
tude larger than the average redshifted 21 cm signal from 
the epoch of reionization. Fortunately, the foreground 
contamination is expected to be spectrally smooth and 



* Since spin temperature fluctuations are expected to be coherent 
on rather large scales |45| . they may in fact make it harder to 
extract the signatures of primordial non-gaussianity compared to 
the simplified case considered here. This should mostly impact 
the early stages of reionization, and so we don't consider the 
effects of spin temperature fluctuations further here. 
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distinguishable from the redshifted 21 em signal which 
should have signifieant frequeney structure (e.g. psf). 
Nevertheless, cleaning this contamination will prohibit 
measuring the signal for long wavelength modes along 
the line of sight, and weaken the ability to measure the 
large-scale power spectrum, precisely where primordial 
non-gaussianity should have its most pronounced effects. 
The impact of foreground cleaning depends on the 
bandwidth over which the contamination is estimated, 
and the precise algorithm used to separate or subtract 
the slowing varying line of sight modes (e.g. [4^ l49j). 
One commonly discussed approach is to fit a low or- 
der polynomial in frequency to each interferometric pixel 
in the Fourier (m, v) plane. The optimal choice for this 
method is to use the lowest order polynomial for which 
the residual foreground power is well beneath the sig- 
nal power. Higher order polynomials are required to fit 
the foregrounds if the fitting is done over larger band- 
widths. Previous work suggests that a cubic polynomial, 
described by iV = 4 coefficients, is adequate for fitting 
foregrounds over a bandwidth of 32 Mhz [501.^ For ref- 
erence, the co-moving length scale corresponding to 32 
Mhz of bandwidth is 



iba„d = 390Mpc//i 



■ B ' 




■0.27" 


1/2 


'1 + z' 
9 


32 Mhz 



1/2 



(21) 



A simple estimate is that a polynomial of order N has 
A^ — 1 nodes, and hence removing a polynomial fit of 
this order removes modes with line-of-sight wavenum- 
ber |A:|| mini < iVTr/Lband 15l|. This suggests that the 
minimum measurable wavenumber, k"^ = k?, -\- k\, after 
foreground cleaning is fcmin = |fc||,minl = A^Tr/Lband- If 
A^ = 4 and B = 32 Mhz, fc^in = 0.032/i Mpc'^ Compar- 
ing with Figure [BJ it appears that adequate foreground 
cleaning will remove the Fourier modes where primor- 
dial non-gaussianity has its largest impact k < 0.02/i 
Mpc~^. This rough estimate illustrates the challenge 
of detecting or constraining primordial non-gaussianity 
with redshifted 21 cm measurements, but more detailed 
calculations are required for a quantitative forecast. In 
particular, there may still be a 'window' of scales where 
the distinctive scale-dependent bias signature can be sep- 
arated from foreground contamination. 



A. Fisher Forecasts 

In order to quantify the prospects for constraining /nl 
in the presence of foreground contamination in more de- 
tail, we use the Fisher matrix formalism. 



1 . No foregrounds 

We first consider the constraints that can be achieved 
in the absence of foregrounds; the effects of foreground 
subtraction will be dealt with below. We calculate con- 
straints in the two-dimensional parameter space p = 
(./nl, b^), where b^ is the gaussian bias. In general, (xi) 
- and additional scale-independent clustering enhance- 
ment terms - might be included as additional parame- 
ters. For simplicity, we instead fix (xi) and stick to this 
two parameter model here, since we expect variations in 
&^ to be more important. The Fisher matrix is given by 



J'ap = ^Tr [C„C-iC^C-i] , (22) 



where C is the covariance matrix of the observed data, 
and a, /3 index the components of p.^° The parameter 
covariance matrix, C^g — (PqP^)i (where the average 
is over many cosmological realizations) can then be ob- 
tained from J^, assuming gaussian statistics, by 



{cpy^ = T, 



(23) 



where we have used the superscript p to distinguish be- 
tween the parameter covariance matrix and the covari- 
ance matrix of the observed data. Constraints on /nl 
and 5^ can then be extracted from C^ in the usual way. 

To compute the covariance matrix of the observed 
data, C, we define a data vector, i5, whose components 
Si represent the observed 21 cm brightness temperature 
fluctuations in the voxels of the experiment. The data 
covariance matrix is related to the observable vector by 
Cij = (SiSj), where again the average is over many cos- 
mological realizations. The covariance matrix Cij can be 
written as a sum of a part due to the signal, C's and a 
part due to noise, Cn. 

To simplify the calculation of the signal covariance ma- 
trix, we ignore any anisotropy of the power spectrum 
(due, for instance, to redshift space distortions), which 
should be a good approximation during most of reioniza- 
tion. We divide the total survey volume into sub-surveys 
with line of sight depth corresponding to a bandwidth of 
30 MHz (as we will discuss in a moment, foreground sub- 
traction will be performed separately in each sub-survey) . 
For simplicity we neglect evolution across the line of sight 
depth of each sub-survey. Although this is an imperfect 
approximation, it should be adequate for our purposes 
since foreground cleaning removes the long wavelength 
modes along the line of sight where evolution should be 
most important. Given these assumptions, the signal co- 



^ This full bandpass may be divided up into smaller chunks, of 
perhaps 6 Mhz width, for power spectrum estimation. This is 
to ensure that the power spectrum evolves minimally across the 
redshift range of a chunk |48|| . 



^^ The fact that the mean of our observable vector is zero allows us 
to write the Fisher matrix in this way. 
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variance matrix in a single sub-survey can be written as 



[C5],, = To{z,)To{z,) (^ 



D{0) 



(27r)3 



(24) 



where we have defined 



In principle, we need to include contributions from 
noise, but we instead consider only Cs here, i.e., we work 
in the cosmic variance limit. Presently, our main aim is 
to explore the fundamental limit imposed by foreground 
cleaning, and so it is appropriate to ignore instrumental 
noise here. In this limit, the normalization factor ro(z) 
drops out, and we then have all the ingredients neces- 
sary to compute J- in the absence of foregrounds. Below 
we consider the effects of foreground subtraction on the 
Fisher matrix. 



&2i(fc,z) = -6^(z) 



3(&g(z) - l)/NLangg^£ 

c^D{z)k^T{k) 



+ {l-{x,)). 



In this equation, 621 (fc, z) is the bias factor that converts 
fluctuations in the underlying density field to fiuctuations 
in a 21-cm map. The form here follows from the term in 
the square brackets of Equation [501 The 21 cm bias is 
negative because large scale overdense regions are pref- 
erentially ionized during reionization - and consequently 
dim in 21 cm - while large scale underdensities remain 
neutral and are hence bright in 21 cm. The To{z) factors 
are dimensionful functions that perform the unit conver- 
sions into 21-cm brightness temperature as in Equation 
[T9I D{z) is the linear growth function, Psp,Sp is the lin- 
ear matter power spectrum at present day, and ipiik) are 
the Fourier transforms of the voxel window functions. 
Here Sb is related to the reionization barrier. The precise 
value of b^{z) and Sb vary with the stage of reionization 
(especially with (x^)) and the reionization model. For 
simplicity, we ignore redshift evolution and use fiducial 
values of b^{z) = 1.5, (xi) ~ 0.3 and Sb — 1-0 which are 
in good agreement with both the analytic and numerical 
results presented above. 

We further assume that our voxels are spherical top 
hats of radius Ar. In detail, 21 cm surveys will have 
higher frequency resolution than angular resolution, but 
primordial non-gaussianity impacts only large scale fiuc- 
tuations, and so the higher frequency resolution will not 
help here. This justifies our use of spherical voxels. This 
leads to 



Mk) 



1 



(4/3)7r(Ar)3 j|j|<a. 
jfc-x. f jijkrAr) 
krAr 



3e 



^3^g-^fe■(^-^,) (26) 
(27) 



where voxel i is located at comoving coordinate .t;, fc^ is 
the radial component of k and ji is the spherical Bessel 
function of the first kind. Substituting into Eq. [M] we 
find 



[CsL = 9To{z,)Toiz,) 



Diz} 
D{0) 



d In k 



_fc3 



■Ps,.5Ak){b2lik,z)y 



jo(k\xi-Xj\) 



Ji(fcAr) 
kAr 



(28) 



2. Effects of Foreground Subtraction 



(25) 

One way to incorporate the effects of foreground sub- 
traction into the Fisher formalism is to add a large 
amount of noise to those modes that arc affected by fore- 
ground removal (i.e. that are most contaminated by fore- 
grounds). See [4^ and [52| for related approaches. Since 
foregrounds are expected to be smooth along the line 
of sight (in the frequency direction) the modes that we 
consider here are low-order polynomials along the line of 
sight and each mode is non-zero in only a single angular 
pixel. In effect, this means that we are assuming the fore- 
grounds are uncorrelated between different angular pixels 
(a conservative assumption). We will define m°'{9i, Zj) to 
be the value of the ath mode in the ith angular pixel and 
the jth redshift bin. So, for instance, the constant mode 
in the angular pixel labeled by a is 



i, Zj) 



1, for all j a i = a, 
0, for all j ii i ^ a. 



(29) 



Higher order modes correspond to higher order polyno- 
mials in Zj. 

To incorporate mode subtraction into the Fisher for- 
malism, we define a constraint matrix 

Ccon{0i,Zj,0i',Zj>j:=^K ^ m°-{0^,Zj)m'^{6i',Zj'), 

a=l 

(30) 

where mf is the value of the ath mode in the ith voxel 
and K is some large number. The data covariancc matrix 
is then adjusted by the constraint matrix: 



C — Cs + Ceo 



(31) 



and the computation of the Fisher matrix proceeds as 
before. 

Ref. [501 -h^^s suggested that foregrounds can be fit with 
a cubic polynomial over a bandwidth of ^ 30 MHz. We 
therefore consider the effects of foreground subtraction 
separately in each subsurvey. We explore the effects of 
foreground subtraction at this level, and also vary the 
number of foreground modes (using higher order polyno- 
mials) to better understand the robustness of constraints 
on /nl to foreground subtraction. 



3. Fisher Results 
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In order to quantify the eonstraints on /nl that ean 
be obtained in the future using the technique presented 
here, we consider the case of a fuh-sky survey that covers 
a frequency range from 120 Mhz to 210 Mhz. This survey 
covers redshifts between z — 5.8 and z — 10.8, centered 
on (z) — 8.3}^ Both the survey area and the large band- 
width of our hypothetical survey are optimistic, but this 
is appropriate for exploring the ultimate limits on /nl 
constraints from the reionization-era 21 cm signal. As 
discussed above, the survey volume is divided along the 
line of sight into three subvolumes, each with depth 30 
Mhz. We divide the survey region into voxels measur- 



ing 



across the sky and 3 Mhz in the frequency 



direction. The resulting voxels are roughly cubical with 
side length (diameter) 50 Mpc so that our spherical voxel 
assumption is not a bad approximation. The effects of 
using finer voxelizations are explored below. 

As our fiducial survey has several million voxels, the 
resulting covariance matrix is very large. Rather than at- 
tempt to invert this large matrix, we compute the Fisher 
matrix for a 5° x 5° region (with equivalently sized vox- 
els) , and scale the resulting Fisher matrix to account for 
greater sky coverage. Computing the Fisher matrix in 
this way assumes that no information is contributed by 
pairs of voxels separated by more than 5°. This is a con- 
servative assumption and we expect it to be a reasonable 
approximation. We assume that foreground removal is 
performed separately for each line of sight and for each 
of the three subvolumes of the full survey. 

Figure [9] shows the projected constraints for the sur- 
vey described above in the two dimensional parame- 
ter space defined hy p = (/nlj^^T)- Fiducial values 
of /nl and 6^ are marked with a cross-hair; they are 
(/nl,&^) = (5,1.5). We show both the constraints ob- 
tained in the absence of foregrounds (red, dotted curves) 
and those obtained when foreground subtraction is im- 
plemented in the manner described above (black, solid 
curve) . Although the foregrounds degrade the detection, 
it seems that significant constraints are still reasonable 
even after foregrounds are subtracted. Specifically, the 
1-(T constraints on /nl after marginalizing over 6^ are 
0.43 and 1.8 before and after removing foreground con- 
taminated modes, respectively. 

To get a better sense of the range of scales that con- 
tribute most to the predicted constraint on /nl in the 
presence of foregrounds, we now consider the effects of 
varying the number of voxels used in the experiment 
and also the number of foreground modes that arc sub- 
tracted. Decreasing the number of voxels effectively de- 
creases krnax , the maximum wavenumber used in the con- 



~r 



~r 



~r 



With foreground mode subtraction 
No foreground mode subtraction 




^^ The Universe is likely fully ionized at the low redshift end of the 
range considered here, but we don't expect our results to depend 
sensitively on the precise redshift range considered. 



FIG. 9: Fisher projections for the error on /nl and b^ with 
and without foreground mode subtraction in the two dimen- 
sional parameter space defined by /nl and be- The cross 
shows our fiducial parameter values. The survey considered 
here is a futuristic, full-sky, cosmic-variance limited 21 cm 
survey covering 120-210 Mhz. 



straint, while increasing the number of foreground modes 
effectively increases kmin- Here, rather than consider an 
experiment across the full bandwidth, we consider one of 
the sub-surveys, ranging from 120 Mhz to 150 Mhz. 

Figure [TOl shows the projected errors on /nl as a func- 
tion of the voxel size of the survey. In generating this 
figure we have assumed the same fiducial values of /nl 
and b^ as above: (/nl, b^) = (5, 1.5). It is clear that the 
Ifj error on /nl declines rapidly with decreasing voxel 
size until the voxel diameter is roughly 50 Mpc; little in- 
formation is gained by using smaller voxels. This means 
that most of the information on /nl is coming from large 
scales, as anticipated. 

We can get a handle on the maximum scale that con- 
tributes to the /nl constraint after foreground cleaning 
by varying the number of foregrounds modes that are 
subtracted. Figure [TT] shows the projected errors on /nl 
as a function of foreground modes that we subtract across 
the 30 Mhz bandwidth. The upper x-axis shows the cor- 
responding maximum scale that can be constrained by 
the data, given by L/N^odcs, where iVmodcs is the num- 
ber of foreground modes subtracted, and L is the dimen- 
sion of the survey. Based on previous studies (e.g. J5C| ) . 
we expect A^ = 4 modes to be sufficient to remove fore- 
grounds over the bandwidth considered. Provided this 
is indeed sufficient, there should be a narrow range of 
scales that arc impacted significantly by /nl, yet sur- 
vive foreground cleaning. In particular. Figures [TU] and 
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With foregrounds, Marginalize bias 
With foregrounds, Fix bias ^ 

No foregrounds. Marginalize bias ^ 
No foregrounds. Fix bias ^ 




65 70 75 80 85 
Voxel Diameter (Mpc) 

FIG. 10: Fisher projections for the error on /nl and 6^ with 
and without foreground mode subtraction as a function of 
the voxel radius of the experiment. The projections here are 
for a survey with 30 Mhz bandwidth rather than the 90 Mhz 
bandwidth assumed in Fig. [9] 



L/N„„,,, (Mpc) 
100 75 



S 6 



Marginalize bias 
Fix bias 



Number of Modes 

FIG. 11: Fisher projections for the la error on /nl as a func- 
tion of the number of foreground modes that are subtracted 
(over a bandwidth of 30 MHz, with 50 Mpc diameter voxels). 
Each foreground mode is a polynomial so that e.g. subtracting 
6 modes corresponds to fitting a 5th order polynomial across 
the bandwidth of the experiment. Based on previous studies 
N = 4 is expected to be sufficient to clean foregrounds. 



[TTl demonstrate that the constraint on /nl comes mostly 
from the narrow range of scales between 50 Mpc to 125 
Mpc, corresponding to roughly k ~ 0.01 — 0.03 hMpc^^ . 
The ability to tightly constrain /nl, in spite of this lim- 
ited dynamic range in scale, is driven by the large volume 
of our futuristic survey. This survey samples many modes 
on the scales of interest and thereby provides small sta- 
tistical errors on the power spectrum. It is encouraging 
that this may, in principle, allow constraints that are 
competitive with - or even slightly better than - Planck 
(as seen in Figure H]). The 1-cr error in Figure IH] is 1.8, 
which compares favorably to the existing 1-cr error from 
Planck of 5.8 [l|- However, because of foreground clean- 
ing, the l/(fc^T(fc)) signature will not be observed over a 
large range in scales and will imprint only a slight excess 
power in the largest scale modes. In effort to ensure that 
this signature is robust to foreground cleaning, one might 
test sensitivity to the number of foreground modes that 
are projected-out, along the lines of Figure [TT] 



VI. DISCUSSION AND CONCLUSIONS 

In this work wc have investigated the effect of primor- 
dial non-gaussianity on the bias of ionized regions dur- 
ing reionization. We have extended the analytic model 
of Furlanetto et. al. Q to the case of non-gaussian initial 
conditions and demonstrated that ionized regions will ex- 
hibit a scale dependent bias. Semi-numeric simulations 
confirm these results. 

We have investigated the constraints that measure- 
ments of the power spectrum of fluctuations in the 21 cm 
radiation from the Epoch of Reionization may place on 
the gaussianity of the initial conditions. We find that fu- 
turistic redshifted 21 cm surveys might allow constraints 
that arc competitive with Planck, in spite of foreground 
cleaning. This is still going to be a challenging endeavor, 
since interesting values of /nl lead to only small changes 
in the power spectrum for the modes that survive fore- 
ground cleaning. Future galaxy surveys will also face se- 
vere systematic challenges in order to robustly measure 
the large scale galaxy power spectrum and constrain /nl 
(e.g. [53 )■ The future prospects also need to be con- 
sidered in the context of the recent Planck constraints, 
which almost close the door on the prospects of detect- 
ing primordial non-gaussianity, and using this to probe 
the physics of inflation. Nevertheless, upcoming galaxy 
surveys and redshifted 21 cm measurements will survey 
tremendous volumes, and can measure the large scale 
power spectrum with high statistical precision. These up- 
coming measurements are interesting in their own right, 
and will proceed regardless of searches for primordial 
non-gaussianity. They still deserve careful scrutiny since 
they should be precise enough to reveal subtle signatures 
from small levels of primordial non-gaussianity. Galaxy 
and 21 cm surveys suffer from different systematic con- 
cerns, and their combination may still provide a powerful 
test of primordial non-gaussianity in the post-Planck era. 
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